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A general formulation is developed to study heat conduction in disordered harmonic chains with arbitrary heat baths that 
satisfy the fluctuation-dissipation theorem. A simple formal expression for the heat current J is obtained, from which its 
asymptotic system-size (N) dependence is extracted. It is shown that the "thermal conductivity" depends not just on the 
system itself but also on the spectral properties of the fluctuation and noise used to model the heat baths. As special cases 
of our heat baths we recover earlier results which reported that for fixed boundaries J ~ 1/N 3/I2 , while for free boundaries 
J ~ 1/AT 1 / 2 . For other choices we find that one can get other power laws including the "Fourier behaviour" J ~ 1/N. 

PACS numbers: 44.10.+i, 05.70.Ln, 05.60.-k, 05.40.-a 



The problem of heat conduction in one-dimensional 
classical systems of interacting particles has attracted a 
lot of attention in recent years (lj. A central issue here 
is determination of the dependence of the heat current J 
on system size N. According to Fourier's law one expects 
J ~ 1 /N but a large number of studies ]l[-p|,p|"^3[ sug- 
gest that in one dimensions this may not always be true. 
Instead one finds that J ~ 1/N a where a is usually differ- 
ent from one. Some obviously interesting and important 
questions are: what are the necessary and sufficient con- 
ditions under which a = 1?, what does a depend on and 
are there universality classes? One of the main problems 
in the field has been that most studies have been lim- 
ited to numerical simulations of nonlinear systems and it 
has been difficult to arrive at definite conclusions from 
results of such studies. Thus, so far no clear understand- 
ing has emerged. We note that numerical simulations are 
problematic because: (i) accurate numerical solutions of 
nonlinear equations are very time consuming, (ii) equili- 
bration times are typically very long and this limits one 
to small system sizes and (iii) dependence on boundary 
conditions not clearly understood. 

One of the earliest models to be investigated was the 
disordered harmonic chain |^ || . This problem is ana- 
lytically tractable to a large extent and the exponent a 
has been obtained analytically, though in a semi-rigorous 
way. Surprisingly a seems to depend on boundary condi- 
tions: for fixed boundary conditions (the Lebowitz model 
j|) a = 1/2, while for free boundaries (the Rubin-Greer 
model ||) a = 3/2. This dependence on boundary con- 
ditions has not been understood in a precise way. 

In this paper we revisit this problem. We present a 
general formulation of the problem which enables one 
to view the two different boundary conditions as two 
special cases of a range of possible thermal reservoirs 
satisfying the fluctuation dissipation theorem. An ap- 
proximate scheme, based on results from the theory of 
product of random matrices, along with inputs from our 
numerical studies, enables us to obtain the asymptotic 



TV-dependence of the current. We find the surprising re- 
sult that the exponent a depends not only on the proper- 
ties of the disordered chain itself, but also on the spectral 
properties of the heat baths. For special choices of baths 
one gets the "Fourier behaviour" a = 1. 

We consider heat conduction through a one- 
dimensional disordered harmonic chain. Particles i = 
1,2...N with random masses are connected by harmonic 
springs with equal spring constants (set to the value 1). 
The Hamiltonian of the system is thus 



H = 



i=i 



yi^ + y 



A' 



(xi - Xi +1 ) 2 
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where {x{\ are the displacements of the particles about 
their equilibrium positions, {pi} their momenta and {mi} 
are the random masses. We put the boundary conditions 
xq = iC/v+i = 0. The particles in the bulk evolve through 
the classical equations of motion while the boundary par- 
ticles, namely particles 1 and N are coupled to heat 
baths. The coupling to heat baths is effected by including 
dissipative and noise terms in the equations of motion of 
the end particles. The choice of the dissipative and fluc- 
tuating forces is not unique. Different forms can be cho- 
sen provided that they satisfy the fluctuation-dissipation 
theorem. 

We consider the following equations of motion for the 
particles: 



m\X\ 
mii'i 
tyinx'n 



where the terms A^^t) and r]T,^(t) describe dissipa- 
tion and noise, and will be specified later. We as- 
sume, unlike B, that the heat baths have been switched 
on at t = — oo. To obtain the particular solution to 



-2 Xl +x 2 + / dt 1 A L (t - t^xxit 1 ) +n L (t) 
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-2x N + xj\[ —i + f dt'A R (t - t')x N (t') + t] R (t) 
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these set of equations we define the Fourier transforms 

= /o°°^i,ii(t)e- 4ttf *. Plugging these into Eq. (§) 
leads to the following particular solution: 



xi(t) = 



1 

2^ 



dw^ m 1 (w)?7 m (w)e 4 ' 



where 



(3) 
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2 tCjT 



Y = $ - urM - A(it>) with 

M im = mi6i im ; A lm = Si im (A L (u>)6is + A R (u)Si tN ) 
fjl = + Vr(u)8i,n- 

The full solution at time t would be the sum of this partic- 
ular solution and a general solution of the homogeneous 
equation, which would depend on the initial conditions. 
Since we are interested in the steady state properties only, 
we will not require the general solution. 

We now specify the properties of the dissipation and 
noise. Let us consider a system driven by a stationary 
noise rj(t) with the following correlator 



(r](cu)f]{uj')) = 2ttTI(uj)S(u + u') 



(4) 



If the dissipation is given by A{uf) = a(uj) — ib(uS), where 
a(u>) and b(uo) are real, then it follows from the fluctua- 
tion dissipation theorem || that the choice 



I{lo) = 2b(u)/w 



(5) 



ensures thermal equilibration of the system to the tem- 
perature T. We choose the same I(oj) and A(u>), satis- 
fying Eq. (H), at both boundaries. The noise correlators 
given by Eq. (Q) are made different by setting T = Tl at 
the left end and T = Tr at the right end. For thermal 
equilibration it is necessary that the range of frequencies, 
over which I{oS) is non-zero, includes the normal modes 
of the disordered chain, and we will only consider cases 
where this is true. 

For any given disorder realization, the energy current 
in the steady state is given by 

J={[[ dtM L (i-0zi(*')+^(*)]zi(*)>, (6) 



where (...) denotes a noise average. Using Eqs. (§,[ 
and after some algebraic manipulations, this reduces to 
the following simple form: 



J 



Tr - T, 



^ /> OO 

— ■ / dujt 2 N {uj) where 

J — oo 

4M = 4&(w)Y{£(w)Y{£(-w) 



(7) 



We note that t 2 N (ui), which is like a transmission coeffi- 
cient, depends both on the system and bath properties. 
We now proceed to write the current in a form where the 
separate effects of the bath and system are more explicit. 
We first note that 



^m^m = lAjvMr 2 with 

A n (lu) =Det[Y] 

D hN - A{uj)(D 2 ^n + £>i,jv-i) + A 2 (uj)D 2 ,n-i 
D\n —Di N-i \ I 1 
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where D^ m is defined to be the determinant of the sub- 
matrix of $ — ui 2 M beginning with the Zth row and col- 
umn and ending with the mth row and column. Clearly 
Di.rn depends on the system alone while A(uj) depends 
on the bath. We further note the following result which 
is easy to prove: 



N 



Ti 



D 2 ,n —D 2 ,n-i 

2 — miuj 2 — 1 
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= T l T 2 ....T 1 
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where 
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The results of j^j] follow from the following choices of 
heat baths: 

(i) Lebowitz model : A(uj) — —i-ycj; I(u) — 2j, (10) 

(ii) Rubin — Greer model : 

AM = 1 - ^ - i|(4 - c 2 ) 1/2 ; J( w ) = (4 - co 2 )'/ 2 \ u \ < 2 



A(c 



1- y + |(4-c 2 ) 1/2 ; ZM=Q H>2 (11) 



Using these in Eq. (0) wc get the heat currents, Jl and 
Jrg, for the two models respectively as: 

/oo 
duJu 2 j N (uj) 
-OO 

j N (cu) = {2 7 2 w 2 

+D 2 hN + jWiDlv^ + Dl N ) + tV^^J- 1 (12) 
Jrg = (4tt)- 1 (T l - T R ) J duu 2 (A - uj 2 )j N (cu) 
j N (u) = {D 2 N + D\ N _ x + + L>2,jv) 2 



+ 2[2(l-w 2 /2) 2 - l]D hN D % 



N-l 



-2(1 - uj 2 /2)(D ltN + D 2 ,jv-i)(I>i,jv-i + D2,n)}-\ (13) 

which are the same as in J|,|| (the differences are due 
to a slightly different convention used by us). Semi- 
rigorous arguments |,|,| indicate that (J L ) ~ l/N 3 / 2 
while (Jrg) ~ where the angular brackets now 

denote a disorder average. For finite chains it is straight- 
forward to numerically compute the integrals appearing 
in Eqs. ( ^2|Jl~^ ) for given realizations of disorder and then 
perform disorder averages to obtain ( Jl) and (Jrg) ■ Wc 
show the results in Fig. (|l|) for the case where the masses 
are chosen from a uniform distribution between 1 — dm to 
1 + 5m. We do get the expected power-law behaviours. 

We now present a scheme which allows us to determine 
the iV-dependence of the current for arbitrary choices of 
heat baths. This is based on the following observations: 
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(i) The first observation follows from the Furstenberg 
theorem on the limiting form of product of random non- 
commuting variables. For the case considered here, the 
theorem states that, for almost any choice of the sequence 
of random masses {mi}, 

lim — log |TiT 2 ....T w u| = 7(w) > (14) 

N—>oc iV 

for any non-zero vector u, with fluctuations of order 
l/y/N. Further it can be shown that |2]j|] in the limit 

7(w) -> (< m 2 > - < m > 2 )w 2 /(8 < m >). (15) 

This means [from Eq. (^)] that the D\ <m , which occur 
in the denominator of the integrand in Eq. (Q), diverge 
exponentially with N , and hence the only significant con- 
tribution to the current comes from low frequency com- 
ponents of order ~ 1/iV 1 / 2 . We note that the fact that 
low frequency modes are extended follows from the trans- 
lational invariance of the random-mass model. 

(ii) The result Eq. (E5|) has been obtained in the strict 
limit of N — > oo when the ratio of successive particle 
displacements reaches a stationary state. For finite N, 
we find from our numerical studies that this result is 
true only for uj ~ 1/N 1 / 2 . In Fig. (||) we have plotted 
(|-Di,tv|) as a function of frequency. We find that the 
exponential growth predicted by Eq. ( |l5| ) does not occur 

atw ~ 1/N 1 ' 2 . In this range we find instead [see Fig. (^)] 
that Di.n is very accurately given by its form for the 
ordered case with masses all equal to < m >= 1. Thus 
over the range u> *~ , we shall approximate j{to) in 

Eq. (^) by its form for the ordered chain. We expect this 
approximation to be good as long as we are interested 
only in the asymptotic iV-dependence. 

For the equal mass case one has -Di,jv = 
sin [k(N + l)]/sin (k) where u> — 2sin(fe/2). Hence 
within our approximate scheme we then get the following 
expression for the disorder- averaged current: 

(J) ~ (T L - T R ) / f(k)dk (16) 
Jo 

f(k) = b 2 {uj) sin 2 (fc) cos(fc/2) x {| sin[fc(7V + 1)] 

- 2A(u) sin(fciV) + A 2 (lu) sin[k(N - 1)]|}~ 2 . 

It is clear that the form of A(oj) at low frequencies will 
determine the asymptotic A^-dependence of the current. 
For the Lebowitz model A(lu) = —ijcu while for the 
Rubin-Greer model A(u>) ~ 1 — iu> and Eq. jl^ ) does give 
the expected 1/N 3 / 2 and 1/N 1 / 2 behaviour for the two 
cases respectively. In general we find that J ~ where 
the exponent a depends on the low-frequency behaviour 
of A(ui). Some special cases are: 

(i) A(u>) ~ —isgn{ui)uj s : Eq. ( pif ) then gives a = s/2 + 1 
for s > 0. 



(ii) A(uS) ~ 1 — isgn(uj)ui s : in this case we get a = 1 — s/2 
for < s < 1 and a — s/2 for s > 1. Note that the 
case s = 2 gives a — 1 that is, a Fourier-like behaviour. 
We verify this by a numerical evaluation of the integral 
in Eq. (^) for chains of finite length and given disorder, 
and then averaging of the current over many disorder 
realizations. The result is shown in Fig. ([!]). 

One simple way of generating thermal sources with 
different spectral properties is to couple the disordered 
chain to an infinite set of oscillators in thermal equi- 
librium. The distribution of oscillator frequencies can 
be arbitrary except that it should include the range of 
the disordered chain frequencies. In this case it can 
be shown that the equations of motion are of the gen- 
eral form Eq. (^) with A(t) — / °° G(uj q ) sin (ui q t)duj q 
where G(u> q ) depends on the choice of oscillator frequen- 
cies. The Rubin-Greer model, where the bath is sim- 
ply an infinite ordered chain, corresponds to the choice 
G(uj q ) = ^w g (4 — to 2 ) 1 / 2 for uj q < 2 and zero elsewhere. 

Finally, we have also studied the effect of introducing 
a quadratic external potential, in addition to the mass 
disorder. In this case, the low frequency current-carrying 
modes are suppressed and we find that the current decays 
exponentially with system-size. 

In summary we have studied the nonequilibrium steady 
state of a mass disordered harmonic chain coupled to heat 
baths at different temperature. We have shown that the 
system size dependence of the energy current, given by 
J ~ 1/N a , is determined not just by the properties of 
the system itself but also by those of the heat baths. 
One gets a continuous set of exponents a depending on 
the low frequency spectral properties of the bath. This 
seems contrary to the general belief in nonequilibrium 
statistical mechanics that the steady state of a close-to- 
equilibrium system will not depend on the details of the 
boundary conditions sustaining the steady state. We ex- 
plain this by arguing as follows: the integrability of the 
harmonic system prevents local thermal equilibrium from 
being established and so the system is really far-from- 
equilibrium. Each of the noninteracting modes indepen- 
dently carries some energy current and the total current 
depends on how exactly the heat baths distribute en- 
ergy among the various modes. For non-integrable sys- 
tems, we expect that there should be transfer of energy 
amongst various modes leading to a state of local thermal 
equilibrium. Hence energy transport should be indepen- 
dent of boundary conditions. However careful studies 
are needed to verify that this is actually true, especially 
for systems which may be close to integrable ones. In- 
deed the possibility of boundary condition dependence 
is suggested from the results of some studies on nonlin- 
ear models SJlQ]. These have often been attributed to 
the fact that some boundary conditions lead to jumps in 
the temperatures at the boundaries which in turn makes 
it difficult to extract the correct system-size dependence 
of current. However our study shows the possibility of 
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boundary-condition dependence even in the absence of 
such jumps. 



Some other interesting questions are: (1) are the pecu- 
liarities of the harmonic chain generic to any integrable 
system, (2) are these results also true for harmonic sys- 
tems in higher dimensions and (3) can the present formal- 
ism be extended to the quantum-mechanical case. The 
answers to these questions may have implications for un- 
derstanding experiments on heat conduction in systems 
such as insulating nanowircs, where similar boundary re- 
lated effects could lead to modification of Landauer-type 
formulas for thermal conductivity 

I thank Madan Rao for very helpful discussions. 
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FIG. 1. System size dependence of the disorder averaged 
steady state current for three different models of heat baths. 
The straight lines shown have slopes 1/2, 1 and 3/2. In all 
cases the disorder strength 8m = 0.2. The error in the mea- 
surements is of the order of the size of the points. 
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FIG. 2. Growth of solutions in a random harmonic chain 
for N = 10 6 and N = 10 4 (inset). The disorder strength 
was taken to be 5m = 0.2. Note that the exponential growth 
starts from oj ~ c/vN (with c « 13). The smooth solid curves 
correspond to the exponential growth predicted by Eq. dl3). 
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